library(plyr)
library(estimatr)
library(emmeans)
library(tidyverse)
library(magrittr)
library(broom)
library(metafor)
library(ggbeeswarm)

d1 <- "https://github.com/thomasjwood/full_fact/raw/main/ff_mt.rds" %>% 
  url %>% 
  gzcon %>% 
  readRDS %>% 
  mutate(
    out_agree = out_agree %>% 
      mapvalues(
        c("Tend to disagree",
          "Tend to agree"),
        c("Disagree",
          "Agree")
      ) %>% 
      factor(
        c("Strongly disagree",
          "Disagree",
          "Neither agree nor disagree",
          "Agree",
          "Strongly agree")
      ),
    out_true = out_true %>% 
      factor(
        c("True",
          "Probably true",
          "Not sure",
          "Probably false",
          "False") %>% 
          rev
      ),
    agree_num = out_agree %>% 
      as.numeric,
    true_num = out_true %>% 
      as.numeric,
    cond = cond %>% 
      factor(
        c("cond_items",
          "cond_misinfo",
          "cond_corr")
      ),
    cond2 = cond
  )


# sequential contrast coding
contrasts(d1$cond2) <-  MASS::contr.sdif(n = 3)

d1$com_num <- d1$agree_num %>% 
  add(
    d1$true_num
  ) %>% 
  divide_by(2)


t1 <- d1 %>% 
  mutate(
    country = issue %>% 
      str_detect(
        "global|saltwater"
      ) %>% 
      ifelse(
        "Multi-country",
        country
      )
  ) %>%
  filter(
    wave == "wave 1"
  ) %>% 
  group_by(
    country, wave, issue
  ) %>% 
  nest 


t1$mods <- t1$data %>% 
  map(
    ~lm_robust(
      com_num ~ cond2, 
      data = .x 
    )
  )

d1 %>% 
  group_by(cond2) %>% 
  summarize(mu = com_num %>% mean)

# cond2           mu
#   <fct>        <dbl>
# 1 cond_items    2.96
# 2 cond_misinfo  3.04
# 3 cond_corr     2.53

t2 <- t1$mods %>% 
  map_dfr(
    ~tidy(.x)
  ) %>% 
  as_tibble %>% 
  filter(
    term %>% 
      str_detect("Intercept") %>% 
      not
  ) %>% 
  mutate(
    type = term %>% 
      str_detect("2-1") %>%
      ifelse(
        "Misinformation effects",
        "Correction effects"
      ) %>% 
      factor(
        c("Misinformation effects",
          "Correction effects") %>% 
          rev
      )
  ) %>% 
  group_by(type) %>% 
  nest

t2$rma <- t2$data %>% 
  map(
    function(i)
      rma.uni(yi = estimate, sei = std.error, data = i)
  )

t2$rma %>% 
  map(plot)

t3 <- t2 %>% 
  pmap_dfr(
    function(type, data, rma)
      
      data %>% 
      mutate(
        type =type
      ) %>% 
      bind_rows(
        rma %>% 
          tidy %>% 
          select(-type) %>% 
          mutate(
            conf.low = estimate %>% subtract(std.error %>% multiply_by(1.96)),
            conf.high = estimate %>% add(std.error %>% multiply_by(1.96)),
            type = type %>% 
              factor(
                c("Misinformation effects",
                  "Correction effects") %>% 
                  rev
                )
            )
        )
    ) %>% 
  mutate(
    sig_insig = p.value %>% 
      is_less_than(.05) %>% 
      ifelse(
        "Significant",
        "Insignificant"
        )
    ) 


t3$lab <- t3$estimate %>% 
  round(2) %>% 
  str_c(
    t3$p.value %>% 
      gtools::stars.pval()
  ) %>% 
  str_squish 

t3$lab %<>% 
  str_extract("\\d") %>% 
  equals("0") %>% 
  ifelse(
    t3$lab %>% 
      str_remove("\\d"),
    t3$lab
  )

t3$lab %<>% 
  str_sub(-1) %>% 
  equals(".") %>% 
  ifelse(
    t3$lab %>% 
      str_sub(, -2),
    t3$lab
  )

t3 %>% 
  ggplot() +
  geom_hline(
    yintercept = 0,
    linetype = 2
  ) +
  geom_dotplot(
    aes(
      type, estimate, fill = sig_insig
    ),
    dotsize = .6,
    stackdir = "center", 
    binaxis = "y", 
    binpositions = "all"
  ) +
  geom_linerange(
    aes(ymin = conf.low, ymax = conf.high, x = type, y = estimate),
    data = t3 %>%
      filter(
        term %>% str_detect("overall")
      ),
    size = .5,
  ) +
  geom_point(
    aes(
      type, estimate
    ),
    fill = "white",
    size = 11,
    shape = 21,
    data = t3 %>%
      filter(
        term %>% str_detect("overall")
      )
    ) +
  geom_text(
    aes(
      type, estimate, label = lab
    ),
    label.size = 0,
    fill = "grey95",
    size = 2.65,
    family = "Roboto",
    data = t3 %>%
      filter(
        term %>% str_detect("overall")
      )
    ) +
  labs(
    y = "Difference on 5pt scale",
    x = "",
    fill = ""
    ) +
  scale_fill_grey(
    start = 1, end = .45,
    breaks = t3$sig_insig %>% 
      unique,
    labels = c("p < .05", "p >= .05")
    ) +
  scale_x_discrete(
    breaks = t3$type %>% 
      levels,
    labels = t3$type %>% 
      levels %>% 
      str_replace(" ", "\n")
    ) +
  coord_flip() +
  theme(
    strip.text.x = element_text(
      size =  11, 
      angle = 0
    ),
    strip.text.y = element_text(
      size =  11, 
      angle = 0
      ),
    axis.text.y = element_text(hjust = .5),
    legend.position = c(.1, .8)
    )
